Principal Component-based Radiative Transfer Model (PCRTM) 
For Hyperspectral Sensors, Part I: Theoretical Concept 


Xu Liu a , William L. Smith b ’ Daniel K. Zhou a , Allen Larar a 
a NASA Langley Research Center, MS401A, Hampton, VA 23681, USA 
b Hampton University, VA 23668, USA 


ABSTRACT 


Modern infrared satellite sensors 1 ' 5 such as AIRS, CrIS, TES, GIFTS and IASI are capable of 
providing high spatial and spectral resolution infrared spectra. To fully exploit the vast amount of 
spectral information from these instruments, super fast radiative transfer models are needed. This 
paper presents a novel radiative transfer model based on principal component analysis. Instead of 
predicting channel radiance or transmittance spectra directly, the Principal Component-based 
Radiative Transfer Model (PCRTM) predicts the Principal Component (PC) scores of these 
quantities. This prediction ability leads to significant savings in computational time. The 
parameterization of the PCRTM model is derived from properties of PC scores and instrument line 
shape functions. The PCRTM is very accurate and flexible. Due to its high speed and compressed 
spectral information format, it has great potential for super fast one-dimensional physical retrievals 
and for Numerical Weather Prediction (NWP) large volume radiance data assimilation applications. 
The model has been successfully developed for the NAST-I and AIRS instruments. The PCRTM 



model performs monochromatic radiative transfer calculations and is able to include multiple 
scattering calculations to account for clouds and aerosols. 
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1. INTRODUCTION 


With the development of modem remote sensing instruments, progressively higher spectral and 
spatial resolution remote sensors are being launched. These modem sensors such as AIRS, TES, 
CrIS, IASI, and GIFTS have thousands of channels and provide a wealth of information on 
atmospheric and surface properties. The new sensors also present challenges to the analysis of 
data. In order to analyze and utilize these vast amounts of data, efficient ways to capture the 
information content of the measurements are needed. One of the key components of analyzing 
these data is the radiative transfer forward model, which relates the atmospheric states (such as 
temperature and trace gas profdes) to the radiances observed from the advanced sensors. Due to 
the high spectral resolution, a large number of Radiative Transfer (RT) calculations through the 
inhomogeneous atmosphere are needed. Therefore, it will be computationally expensive to 
incorporate these data into global NWP data assimilation models using conventional radiative 
transfer models. Oftentimes, only subsets of channels are used with existing RT models in order 
to handle the corresponding computational constraints. 


It is well known that performing line-by-line radiative transfer calculations through the 
inhomogeneous atmosphere is a very time consuming process. Most molecular gases in the 
atmosphere have numerous vibrational-rotational transitions, or pure rotational transitions, in the 
infrared spectral region. Molecular line intensities and shapes are non-linear functions of vertical 
profiles of atmospheric temperature, pressure and trace gases. Typically one has to divide the 
inhomogeneous atmosphere into numerous thin layers and use properly weighted atmospheric 
properties within these layers for the radiative transfer calculations. Since the line widths are 
very small in the upper atmosphere, the monochromatic calculations are typically done in 0.0001 to 
0.0008 wavenumber steps in the infrared spectral region. For example, the TES instrument 2 ’ 5 ' 9 
spectral coverage ranges from 650 cm' 1 to 3050 cm' 1 with an average frequency grid of 0.0004 cm' 1 ; 
thus, about 6 million monochromatic radiative transfer calculations are needed to model the whole 
spectral range. 

There are several ways to minimize the computational time needed to perform radiative transfer 
calculations for the processing of hyperspectral data. One of the approaches, which the TES 
science team uses, is to store absorption coefficients as a function of atmospheric pressure and 
temperature at a monochromatic frequency grid in a large lookup table. As a result, the optical 
depths in a particular atmospheric layer can be calculated simply by interpolations and additions. 
In this approach, time-consuming calculations of spectral line shapes and intensities are avoided, 
but the forward model still has to perform numerous monochromatic radiative transfer calculations 
through inhomogeneous atmospheres in order to obtain Top Of Atmosphere (TOA) radiances. To 
minimize the calculations needed, the TES science algorithm takes advantage of the high spectral 



resolution and high information content of the TES measurements. It selects narrow micro- 
windows to perform retrievals for a specific trace gas. However, the approach may be sub-optimal 
in that more channels cannot be used in a simultaneous solution for all desired geophysical 
parameters. 


Another forward model approach is to predict effective layer optical depths by using an efficient 
fast parameterization. The effective layer optical depth is a channel-averaged quantity, which 
contains the Instrument Line Shape (ILS) function or Sensor Response Function (SRF); therefore 
only one radiative transfer calculation per channel is needed. The effective optical depth is 
derived in such a way that the additive property of optical depth (or the multiplicative properties of 
transmittances) between different atmospheric layers and different atmospheric gases holds. 
Equation 1 is an example of how the layer optical depth is calculated: 
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Where x is the effective channel layer optical depth, (j) is a normalized SRF, Av is the spectral span 
of the SRF, / is the atmospheric layer index, t Mono and t ch are monochromatic and channel 
transmittance (space to layer), respectively. T z is pressure weighted layer temperature ratio above 
layer /, and T r is the ratio of the layer average temperature to a layer reference temperature. 6 is 
the satellite zenith angle and P is the atmospheric pressure. 



There are several fast model parameterizations based on effective optical depth for satellite remote 
sensing applications ’ ' . Optical Path Transmittance (Optran), Stand-Alone Radiative Transfer 
Algorithm (SARTA) and RTTOV are three well-known fast models of this kind 26 27 37 . These 
approaches have been successfully used in the operational algorithm for HIRS, AMSU, and AIRS 
data retrievals. 

Another method of fast parameterization is to predict channel transmittances or radiances by using a 
few representative monochromatic transmittances or radiances l0 " l6 ‘ l9 ‘ 2l ‘ 24, 3 1 ‘ 34 ‘ 35 . Correlated K 
Distribution (CKD), Exponential Sum Fitting Transmittance (ESFT), Radiance Sampling Method 
(RSM) and Optimal Spectral Sampling (OSS) are examples this type of fast RT model approach. 
Equations 2 and 3 show how the channel radiances or transmittances are calculated. 
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In the equations, R cl, (v) and l <h (v) represent channel radiance and transmittance at a center 
frequency v , respectively, and i-v, is the weight for the pre- determined monochromatic radiance or 
transmittance. Both CKD and ESFT have the advantage of being monochromatic (can extend to 
include multiple scatterings) and very fast (only a few monochromatic calculations are needed per 
channel). They are also very accurate for a single atmospheric layer and a single absorption gas. 



However, these models are usually trained on one atmospheric layer, and the dependency of 
channel transmittance on pressure, temperature and gas amounts are introduced later on by 
assuming a good correlation between vertical layers and no correlation between overlapping gases. 
Extending them to vertically inhomogeneous atmospheres leads to limited accuracy. Extensive 
efforts were made in the scientific community to extend these methods to remote sensing 
applications 1M3 ' 15 ' 24 . One way to overcome this deficiency is to include both layer transmittance 
and space-to-layer transmittance in the training 12 , and another way is to predict the TOA channel 
radiance using statistically selected TOA monochromatic radiances 11 (i.e. RMS method). This has 
the advantage of treating overlapping gases and inhomogeneous layers consistently for all 
atmospheric layers. The drawback of the RSM is that it needs a relatively large number of 
monochromatic radiances to approximate TOA channel radiances due to the method used to 
determine the monochromatic frequency locations. The OSS method developed by AER 34 " 36 
overcomes this deficiency by fitting TOA channel radiances using a robust ESFT minimization 
scheme 10 , where typically, one to fifteen monochromatic radiative transfer calculations are needed 
to predict channel radiance. OSS is also capable of predicting layer or space-to-layer 
transmittances. Moncet et al 34 35 have extended the search method to a more robust “Monte 
Carlo” approach. The OSS fast model is planned to be used to process Cross-track Infrared 
Sounder (CrIS) and Advanced Technology Microwave Sounder (ATMS) data from future 
NPP/NPOESS satellites (http://npoesslib.ipo.noaa.gov/atbd/ cris_atbd_03_09_01.pdf). 

All of the fast models described above are channel-based forward models. These models predict 
channel effective layer optical depths, channel transmittances or channel radiances either by using 
non-linear functions of atmospheric temperature and gas profiles or through reduced 



monochromatic calculations. The PCRTM, on the other hand, does not predict channel radiance 
directly. Instead, it predicts PC scores, which have much smaller dimension as compare to the 
number of channels 14 ’ 38 ’ 39 . The PCRTM will take advantage of the orthogonal properties of the 
principal components (PCs) to compress the spectral information (e.g. channel radiance, 
transmittance or reflectance) into relatively few significant PCs as compared to the number of 
channels. The observed TOA radiances can be converted to PC scores by projecting the spectrum 
onto a set of PCs. Conversely, channel radiances can be calculated by linearly combing PCs using 
PC scores. 


2. THEORETICAL BASIS OF PCRTM FORWARD MODEL 


For high-resolution infrared spectra, there are a lot of channels that have similar properties. For 
example, ozone has more than 100000 absorption lines in the spectral range from 600-3000 cm" 1 . 
Many of these lines have similar Lorentz or Doppler half-widths. Dependencies of the line 
intensity and half-width on atmospheric temperature and pressure are similar. All monochromatic 
radiances are a nonlinear function of atmospheric temperature, moisture and trace gas profiles. 
The number of independent pieces of information is much less than the number of monochromatic 
radiances, therefore many radiative transfer calculations are redundant. The PCRTM approach 
removes this redundancy. The PCRTM will compress channel transmittances or radiances into a 
set of orthogonal eigenvectors or PCs in a descending order according to the associated eigenvalues. 
These radiative transfer properties can be easily regenerated via a linear combination of pre- 
determined significant PCs. The linear coefficients are called principal component scores. These 
PC scores can be obtained by projecting the channel radiances onto each of the PCs. The number of 



significant PCs is determined by checking the RMS errors between the original radiances and 
radiances regenerated using various numbers of PCs. Usually 100-250 PCs are enough to 
regenerate channel radiances to 0.01 K accuracy for hyperspectral sensors with thousands of 
channels. These significant PCs capture spectral variations of TOA radiances from one channel to 
another, while the corresponding PC scores capture the dependencies of spectra on atmospheric 
temperature and gas profiles. Instead of performing monochromatic calculations at thousands or 
millions of frequencies in real time, these RT calculations were performed off-line in a training 
process. Just like other fast forward model trainings 27 ’ 37 , an ensemble of atmospheric profiles, 
which span the expected range of variability, were selected from a large dataset such as TIGR 
profile or ECMWF analyses model outputs 43 . Optical depths at each atmospheric layer were 
calculated at a monochromatic grid of 0.0001 cm' 1 . Radiance spectra were generated using layer 
optical depths under various observation geometries and surface conditions. These 
monochromatic radiances were then convolved with appropriate sensor response or instrument line 
shape functions to form channel radiances. A data matrix whose columns consist of N spectra and 
whose rows consist of channel radiances at M channel frequencies was then formed and a Singular 
Value Decomposition (SVD) was performed on this matrix. The first PC generated by SVD 
represents an average channel spectrum. The second PC is orthogonal to the first one and accounts 
for a major fraction of the variance in the data matrix. Each successive PC is responsible for 
smaller fraction of the total variance in the data. As mentioned earlier, only a few hundred PCs 
are needed to regenerate any spectrum in the data matrix. Since these spectra were generated 
using representative profiles spanning the expected variability of atmospheric and surface 
conditions, the PCs should be able to represent radiances corresponding to other general 
atmospheric and surface conditions. Once the PCs are generated, they are pre-stored in the 



forward model. The task for the off-line training is to find a good way to predict PC scores, which 
in combination with PCs produce channel radiances. For a given atmospheric state, PC scores can 
be predicted using a non-linear function of atmospheric state such as layer average temperature, 
pressure, water layer amounts etc. As derived below, a better way to predict PC scores is to use a 
linear function of a few monochromatic transmittances or radiances. It has the advantage that the 
radiative transfer equation transforms the atmospheric temperature, moisture and trace gas profiles 
into radiances and these radiances have a linear relationship with the PC scores: 
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In this equation, R ch is the channel spectrum vector, U i is the i th PC vector, N pc is the number of 
significant PCs. Y, is the PC score, which is generated by linearly combining monochromatic 
radiances and aj are the associated weights. £ is a vector containing forward model errors. 

Usually, the values of £ are smaller than 0.04 K. The PC score (7 ; ) is a dot product of PC vectors 
(U T Nchxi) with the channel radiance vector ( R ch Nchxi )• 
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The index i represents the ith PC score, j represents the channel index, the superscript T stands for 
transpose, and u(j,i) is an element of the PC matrix (U). The channel radiance is a linear 
combination of monochromatic radiances within the frequency span of that channel. The weights 



are simply the normalized ILS or SRF at the monochromatic frequency grid. Both the pre- 
calculated PC vectors and the ILS (or SRF) do not vary from one atmospheric state to another, 
therefore, the PC score is a linear function of monochromatic radiances (equation 6). 
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where R mono represents monochromatic radiance, </> represent normalized ILS (or SFR) and a is the 
weight. As mentioned before, there is a lot of redundant information among thousands to millions of 
monochromatic radiances or transmittances. There is no need to use all of them to predict 7, in 
equation 6. Unlike methods described above, the PCRTM method selects the location of 
monochromic frequencies in a very straightforward way by using a correlation function 36 . Instead 
of performing a selection within each channel, as is done within OSS, and ESFT models, the 
PCRTM selection is done in the whole spectral range, thus removing maximum redundant 
information in the process. The correlation function correctly identifies redundant information; 
therefore there is no need to perform a labor-intensive search within millions of possible locations. 
The upper panel in Figure 1 shows an example of correlation coefficients as a function of 
frequency. The correlation coefficients are converted to vector angles via an arccosine function. 
The bottom panel is a plot resulting from re-arranging the correlation coefficients according to their 
magnitudes. The monochromatic radiances were selected by choosing predictors with equal 
distance in the values of the correlation coefficients. If the correlation coefficients are equal 
between two monochromatic frequencies, then the corresponding radiance vectors are parallel to 
each other and contain the same information content. Only one of the points is needed. On the 



other hand, if two radiance vectors are orthogonal to each other (i.e. with arccosine of correlation 
coefficient equals to 90 degree), both should be selected as predictors. The weights (a/) of the 
selected predictors (monochromatic radiances) are obtained by a simple SVD regression procedure. 
Usually, only 300-900 terms are needed to predict Y, (i.e. for an infrared spectrum ranging from 
650 to 3000 cm" 1 ). The correlation function method in some way is similar to the CKD method in 
the sense that they both group monochromatic frequencies with similar properties together. In 
the CKD method, monochromatic frequencies inside the spectral span of an ILS or SRF are re- 
ordered according to values of molecular absorption coefficients; while in the correlation function 
method, monochromatic frequencies are re-ordered according to values of correlation functions and 
it is done over the whole instrument spectral range. The increase in the speed of the PCRTM 
model comes from two factors: (1) there are less predictants ( Y,) to predict due to the PC 

compression, and (2) there are less overall predictors due to the use of correlation functions in 
selecting monochromatic frequencies. If the retrieval system uses only PC scores (a few hundreds) 
not channel radiances (a few thousands), the first factor can lead to an order of magnitude saving in 
computational time. The second factor is very significant since the PCRTM only performs 
radiative transfer calculations at a few hundred monochromatic frequencies, while the channel 
based RTMs perform RT calculations at a few thousand frequencies. 

An advantage of the PCRTM model is that no additional effort is required to model spectra with 
different apodization functions. Since the ILS of a weakly apodized spectrum has ringing and 
wide ILS spectral span for each channel, the channel based fast RT models 32 " 34 are more difficult to 
handle as compared to localized ILS. For example, a sine ILS (or boxcar apodization) will 
produce negative channel transmittances. One limitation of an OPTRAN type of model is that a 


logarithm of a negative effective transmittance (see equation 1) cannot be calculated. For models 
using monochromic radiances or transmittances as predictors, you have to choose more predictors 
in order to cover the wide ILS span. This will lead to increased computational time. For 
PCRTM, the instrument ILS information is captured by the pre-stored PCs. Selection of the 
predictors is based on the whole instrument spectral coverage to start with, no extra predictors are 
needed for a weakly apodized spectrum as compared to a strongly apodized spectrum. Another 
advantage of the PCRTM is that it may be practical to include multiple scattering calculations since 
a minimum number of monochromatic radiative transfer calculations are needed. The bulk of the 
multiple scattering calculations are done off-line in the training process. It should be mentioned 
that for the purpose of retrieval, the PC scores contain all of the essential information on the spectra. 
The inversion of state vectors can be done in PC domain without having to go to spectral domain. 
The capability of going to spectral domain is provided mainly for the quality control and for more 
convenient physical interpretation of the retrieval results. 


3. DESCRIPTION OF PCRTM AND PRELIMINARY RESULTS 


A flow diagram of the PCRTM model is shown in Figure 2. The model input parameters include 
eigenvectors for channels spectra, regression coefficients, and absorption coefficients for various 
atmospheric gases. Since the radiative transfer calculations used to generate predictors for PC 
scores are done monochromatically, the radiative transfer calculation is very simple and 
straightforward. In addition to calculating PC scores and channel radiances, the PCRTM is also 



capable of providing an analytical Jacobian matrix, which can be used for one-dimensional (1-D) 
physical retrievals or for 3-D or 4-D variational data assimilations. 


The following example shows how the upwelling radiance can be calculated in a clear atmosphere. 
First, the upwelling radiance ( R " p ) from surface is initiated: 

K P = e v B v (T s ) (7) 

In this equation, £ v is surface emissivity at frequency v, and B v is the Planck function at frequency 
v and at a given temperature. In this case, the surface thermal emission is calculated at surface 
skin temperature T s . The upwelling radiance component at sensor altitude and its derivatives with 

respect to layer temperatures ( v ) and trace gas layer amounts (e.g. v ) can be calculated 
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Here, l represents the layer index, and t, and t nT j_ ; represent layer transmittance and 

accumulative top layer to current layer transmittance, respectively; 6 is the sensor zenith angle. 
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and —1 are the derivatives of layer optical depth with respect to layer temperature and 
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layer water column amount, respectively. — is the derivative of upwelling radiance with respect 

dr, 

to layer optical depth. The superscript 0 indicates that the transmittance is calculated at a zero zenith 
angle (i.e. at nadir direction). These two quantities can be easily computed using a lookup table 
approach 5 . Once the monochromatic radiance and associated derivatives are calculated, equation 
6 is used to predict the PC scores. If channel radiances are needed, equation 4 can be used to 
transform information from PC domain to radiance domain. 


The PCRTM forward model has been implemented for both AIRS and the NPOESS Airborne 
Sounder Testbed Interferometer (NAST-I) instruments. The PCRTM fast model has very high 
accuracy as compared to the line-by-line calculations. Figure 3 shows the RMS and bias errors 
between the LBLRTM 22 and the PCRTM for AIRS instrument for 52 ECMWF diverse profiles at 
different satellite observation angles. 120 PCs and 246 predictors were used for the AIRS PCRTM 
model. The RMS errors for the model are usually less than 0.04K. 

Figure 4 shows the same model comparison for the NAST-I instrument 40 ' 41 . There are 8632 
channels for the NAST-I instrument; only 200 PCs are needed to represent these channels. To 
predict the associated 200 PC scores, only 305 monochromatic RT calculations are needed. This 
means that an average of 0.04 predictors are needed for PCRTM to predict each NAST-I channel. 



In contrast, at least one RT calculation and one predictor are needed for each channel for channel- 
based fast models, requiring a minimum of 8632 RT calculations and predictors. For OPTRAN or 
SARTA type of fast RTM, each channel needs 10 to 30 predictors to predict the effective channel 
transmittances. For ESFT, CKD or OSS type fast RTM, 1 to 15 predictors are needed to predict 
each channel radiance. The PCRTM should have a significant speed advantage as compared to 
channel-based RTMs. 

To test the applicability and accuracy of the PCRTM model, we have compared the observed AIRS 
radiance with radiance calculated using the PCRTM model. The AIRS spectrum was taken from a 
clear sky overpass of Aqua over the Atmospheric Radiation Measurement Tropical Western Pacific 
(ARM-TWP) site on 8 December 2002. The surface pressure and ozone profile are from the 
European Center for Medium-Range Weather Forecasts (ECMWF) model. Surface temperature is 
from University of Maryland at Baltimore County (UMBC) retrieval. Water and temperature 
profiles are mainly from collocated radiosondes over the ARMS-TWP site. The water profile 
above 200 mb and the temperature profile above 60 mb are taken from the ECMWF model 
analyses. The surface emissivity spectrum is obtained from the JHU 42 database. Figure 5 shows 
the observed AIRS radiance (top panel), the PCRTM calculated radiance (middle panel) and the 
difference between the two (bottom panel). Large residuals exist near 1050 cm' 1 . It is clear that 
the ozone profile from ECMWF analysis does not represent the true atmospheric ozone state. The 
spikes on the top and bottom panels are due to the undetected “popping” noise of the AIRS 
instrument detectors. Overall, the agreement is good between the PCRTM calculation and the 
AIRS observation. It should be mentioned that for this version of the PCRTM, the AIRS SRF is 



not the latest available. The channeling and lastest changes in the SRF are now being included in 
the AIRS PCRTM. 


4. DISCUSSION AND CONCLUSIONS 


PCRTM is much faster than channel-based RTMs since it predicts PC scores (or “Super Channel” 
magnitudes) directly instead of channel radiance or transmittance individually. The choice of 
using radiances or transmittance as predictors is derived theoretically from the properties of PC 
scores, where the parameterization of PCRTM is physical-based. The accuracy of PCRTM 
relative to LBL codes is very good, typically less than 0.04K. The accuracy can be improved 
monotonically to reach that of LBL calculations by increasing the number of predictors. Based on 
the application of PCRTM to the AIRS and NAST-I instruments, the PC scores can be predicted 
accurately using a few hundred monochromatic radiances. Channel radiances can be calculated by 
linearly combining pre-stored PCs with PC scores as weights. The radiance variation as a function 
of T, H 2 0, O 3 , CH 4 , N 2 0, CO, surface emissivity and observation geometry is captured via 
monochromatic radiative transfer calculations. The redundant spectral information is captured via 
PC representation. Since RT calculations are done monochromatically, Jacobians can be 
generated with little extra effort and multiple scattering calculations can be included. The PCRTM 
provides both PC scores and the associated Jacobian, therefore it is recommended that the physical 
inversion of state vector be done in the PC domain directly. Information from all channels is 
transformed into PC scores and the retrieval process can take advantage of maximum information 
content and retain the best signal-to-noise ratio of the observed spectrum. This way, there is no need 
to select a sub-set of channels for the sake of computational speed limitation. For quality control, 



channel radiances can be generated with a simple PC transformation. Due to its fast speed and high 
accuracy, PCRTM has great potential for the assimilation of the complete spectrum of observed 
radiances in NWP. Furthermore, the use of the PCRTM may enable cloud parameters in the NWP 
operation. 

This is one of a series of papers we plan to write for the PCRTM model. The main aim of this 
paper is to present the theoretical basis of the method and provide some preliminary results. A more 
detailed description of the model and more results will be presented in future publications. 
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Figure 1. Flpper panel: Arcccosine of correlation coefficients as a function of frequency, bottom 
panel: re-arranged according to the magnitudes. 

Figure 2. Flow diagram of PCRTM 

Figure 3. RMS errors and bias errors of the AIRS PCRTM model relative to FBFRTM 
Figure 4. RMS errors and bias errors of the NAST-I PCRTM model relative to FBFRTM 

Figure 5. Comparison of the observed AIRS radiance over ARM site on 12/08/2002 with radiance 
calculated using PCRTM model. 
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